* Set file directory path using <cd>

* Moran's I is calculated using Maurizio Pisati's spatcorr command, which needs to be installed
net install sg162

********************************************************************************
* To get distance band in kilometer, we used the UTM30N projection, 
// where x and y coordinates are expressed in kilometers (x_km y_km).
// this is for spatcorr, where "x-coordinate and the y-coordinate must be expressed in projected units, e.g., meters, kilometers, miles, or arbitrary digitizing units."

set matsize 11000

use "Appendix Figure A2\data_cross_section_ghana.dta", replace
merge m:1 latitude longitude using "Appendix Figure A2\Ghanagrid_UTM30N"
drop _merge
tempfile moran
save `moran', replace


********************************************************************************
// Spatial matrix needs to be based on the cross-section
// Estimate for each cross-section separately: 1850, 1875, 1897, 1932

**************** Without lagged DV
* 1850 
use `moran', replace
drop if year!=1850 

foreach X in 50 {
foreach Y in 0 {

local geo "prec_mean alt_mean alt_sd soil5 histo_malaria"
local ecogeo "port1850_yn ldist2coast nongcc1873 dist2riv_navigable10 traderoute18501890_dist10 exploroute188089_dist10 dist2rail193210 dist2class123_193010"
local pop "city1800 map08_yn national headchief lupop_1891 lupop_1901 lrpop_1901 lupop_1931 lrpop_1931"
local econ "ln_export_area slavemarket_strict_dist`X' palmoil1900_dist`X' kolanut`Y' rubber1900_dist`X' cocoa minet_dist1932_`X'"

acreg missions_yn `geo' `ecogeo' `pop' `econ' area_sqkm if year==1850, latitude(latitude) longitude(longitude) dist(100) spatial
}
}

predict predicted
gen residual1850= missions_yn-predicted
spatcorr residual1850 , bands(0(15)105) xcoord(x_km) ycoord(y_km) graph
putexcel set "Appendix Figure A2\moran", replace
putexcel A1 = "residual1850" 
putexcel B1 = matrix(r(Moran)), colnames


**************** With lagged DV
* 1875 
use `moran', replace
drop if year!=1875

foreach X in 50 {
foreach Y in 0 {

local geo "prec_mean alt_mean alt_sd soil5 histo_malaria"
local ecogeo "port1850_yn ldist2coast nongcc1873 dist2riv_navigable10 traderoute18501890_dist10 exploroute188089_dist10 dist2rail193210 dist2class123_193010"
local pop "city1800 map08_yn national headchief lupop_1891 lupop_1901 lrpop_1901 lupop_1931 lrpop_1931"
local econ "ln_export_area slavemarket_strict_dist`X' palmoil1900_dist`X' kolanut`Y' rubber1900_dist`X' cocoa minet_dist1932_`X'"

acreg missions_yn missions_yn1850 `geo' `ecogeo' `pop' `econ' area_sqkm  if year == 1875, latitude(latitude) longitude(longitude) dist(100) spatial
}
}

predict predicted
gen residual1875= missions_yn-predicted
spatcorr residual1875, bands(0(15)105) xcoord(x_km) ycoord(y_km) graph

putexcel A9 = "residual1875"
putexcel B9 = matrix(r(Moran))


* 1897
use `moran', replace
drop if year!=1897

foreach X in 50 {
foreach Y in 0 {

local geo "prec_mean alt_mean alt_sd soil5 histo_malaria"
local ecogeo "port1850_yn ldist2coast nongcc1873 dist2riv_navigable10 traderoute18501890_dist10 exploroute188089_dist10 dist2rail193210 dist2class123_193010"
local pop "city1800 map08_yn national headchief lupop_1891 lupop_1901 lrpop_1901 lupop_1931 lrpop_1931"
local econ "ln_export_area slavemarket_strict_dist`X' palmoil1900_dist`X' kolanut`Y' rubber1900_dist`X' cocoa minet_dist1932_`X'"

acreg missions_yn missions_yn1875 `geo' `ecogeo' `pop' `econ' area_sqkm  if year == 1897, latitude(latitude) longitude(longitude) dist(100) spatial
}
}

predict predicted
gen residual1897= missions_yn-predicted
spatcorr residual1897 , bands(0(15)105) xcoord(x_km) ycoord(y_km) graph

putexcel A16 = "residual1897" 
putexcel B16 = matrix(r(Moran))


* 1932
use `moran', replace
drop if year!=1932

foreach X in 50 {
foreach Y in 0 {

local geo "prec_mean alt_mean alt_sd soil5 histo_malaria"
local ecogeo "port1850_yn ldist2coast nongcc1873 dist2riv_navigable10 traderoute18501890_dist10 exploroute188089_dist10 dist2rail193210 dist2class123_193010"
local pop "city1800 map08_yn national headchief lupop_1891 lupop_1901 lrpop_1901 lupop_1931 lrpop_1931"
local econ "ln_export_area slavemarket_strict_dist`X' palmoil1900_dist`X' kolanut`Y' rubber1900_dist`X' cocoa minet_dist1932_`X'"

acreg missions_yn missions_yn1897 `geo' `ecogeo' `pop' `econ' area_sqkm  if year == 1932, latitude(latitude) longitude(longitude) dist(100) spatial
}
}

predict predicted
gen residual1932= missions_yn-predicted
spatcorr residual1932 , bands(0(15)105) xcoord(x_km) ycoord(y_km) graph
putexcel A23 = "residual1932" 
putexcel B23 = matrix(r(Moran))




******************************************************************************
** Figure Spatial correlogram of Moran's I 
import excel "Appendix Figure A2\moran.xlsx", sheet("Sheet1") firstrow case(lower) clear
rename residual1850 residual
replace residual="1850" if _n==1
replace residual=substr(residual,-4,.)
destring residual, replace force
replace residual=residual[_n-1] if  residual==.

* Calculate CIs
gen ci_upper=stat+1.96*sd 
gen ci_lower=stat-1.96*sd

** To offset the CIs
forval x=-0(2)6 {
gen upper_`x'=upper+`x'-3
}

* Figure of Moran's I with lines
twoway (connected stat upper_0 if residual==1850, lstyle(solid) mcol(gs10) lcol(gs10)) (rcap ci_upper ci_lower upper_0 if residual==1850, lcol(gs10)) ///
 (connected stat upper_2 if residual==1875, lstyle(solid) mcol(gs7) lcol(gs7)) (rcap ci_upper ci_lower upper_2 if residual==1875, lcol(gs7)) ///
 (connected stat upper_4 if residual==1897, lstyle(solid) mcol(gs2) lcol(gs2) msymbol(X)) (rcap ci_upper ci_lower upper_4 if residual==1897, lcol(gs2)) ///
 (connected stat upper_6 if residual==1932, lstyle(solid) mcol(black) lcol(black)) ///
 (rcap ci_upper ci_lower upper_6 if residual==1932, lcol(black) ), ///
  yline(0) ytitle("Moran's I") xtitle("Distance bands (km)") ///
 xlabel(15 "0-15" 30 "15-30" 45 "30-45" 60 "45-60" 75 "60-75" 90 "75-90" 105 "90-105")  legend(order(1 "1850" 3 "1875" 5 "1897" 7 "1932") col(4) on)  
graph export "Appendix Figure A2\Appendix Figure A2.pdf"
